home *** CD-ROM | disk | FTP | other *** search
/ MacFormat UK 55 / MF_UK_55_1.iso / mac / Shareware Plus / Developers / VideoToolbox / (Utilities) / MeasureMTF / MeasureMTF.c next >
Encoding:
C/C++ Source or Header  |  1997-07-02  |  11.3 KB  |  357 lines  |  [TEXT/CWIE]
  1. /*
  2. MeasureMTF.c
  3. Copyright © 1990-1993 Denis G. Pelli
  4. Measures the modulation transfer function of a monitor. Uses drifting gratings,
  5. both horizontal and vertical, of constant temporal frequency. The results
  6. are saved in a KaleidaGraph text file, MTF?.data, where ? stands for the screen
  7. number. The data file includes the dc and second harmonic.
  8. HISTORY:
  9. 10/11/90    dgp    wrote it.
  10. 10/12/90    dgp    added dc and harmonics.
  11. 3/18/91        dgp    updated to work with new PlotXY, but not tested
  12. 4/15/91        dgp Check for NewPaletteManager().
  13. 8/24/91        dgp    Made compatible with THINK C 5.0.
  14. 8/27/92    dgp    replace SysEnvirons() by Gestalt()
  15. 10/23/92 dgp try to read latest LuminanceRecord
  16. 2/7/93    dgp    updated to use SetPixelsQuickly.
  17. 7/1/97    dgp    After GDOpenWindow, now call ShowWindow, since GDOpenWindow no longer does.
  18. */
  19. #include "VideoToolbox.h"
  20. #include "Luminance.h"
  21. //#include <Packages.h>
  22. #if (THINK_C || THINK_CPLUS || SYMANTEC_C)
  23.     #include <profile.h>    /* for timing */
  24.     #define SYMANTEC_C_PROFILE 1
  25. #endif
  26. #define TWO_PI (2.0*PI)
  27. #define MAX(a,b) ((a)>(b)?(a):(b))
  28. #define MIN(a,b) ((a)<(b)?(a):(b))
  29.  
  30. void MeasureMTF(void);
  31. double MeasureContrast(GDHandle device,LuminanceRecord *LP,int tPeriod,double c[10],WindowPtr plotWindow);
  32. double saw(double x);
  33.  
  34. void main(void)
  35. {
  36.     assert(StackSpace()>1000);
  37.     StackGrow(30000);
  38.     Require(gestalt8BitQD);
  39.     MeasureMTF();
  40. }
  41.  
  42. void MeasureMTF(void)
  43. {
  44.     GDHandle device,oldGDHandle;
  45.     CWindowPtr window;
  46.     WindowPtr plotWindow,mtfWindow;
  47.     LuminanceRecord LR;
  48.     char string[100],string2[100],outName[100];
  49.     unsigned long seconds;
  50.     FILE *outfile[2];
  51.     int i,tPeriod,width;
  52.     Rect r,srcRect,dstRect;
  53.     double a,b,f,fNominal,c,L,cOut,cH[10],cV[10];
  54.     PlotXYStyle mtfStyle[2];
  55.     unsigned long row[1048];    // hopefully long enough for any video device
  56.     short oldScreen;
  57.  
  58.     assert(StackSpace()>4000);
  59.     MaximizeConsoleHeight();
  60.     #if (THINK_C || THINK_CPLUS || SYMANTEC_C)
  61.         console_options.ncols=100;
  62.         printf("\n");                       /* Initialize QuickDraw */
  63.     #else
  64.         InitGraf(&qd.thePort);
  65.         InitFonts();
  66.         InitWindows();
  67.         InitCursor();
  68.     #endif
  69.     #if SYMANTEC_C_PROFILE
  70.         InitProfile(200,1);
  71.         _profile=0;    /* disable profiling for the moment */
  72.     #endif
  73.     printf("Welcome to MeasureMTF.\n");
  74.     #include "LuminanceRecord1.h"                    // read at compile time
  75.     oldScreen=LR.screen;
  76.     if(GetScreenDevice(1)!=NULL){
  77.         printf("Which screen would you like to calibrate (%d):",LR.screen);
  78.         gets(string);
  79.         sscanf(string,"%d",&LR.screen);
  80.     }
  81.     else LR.screen=0;
  82.     sprintf(string,"LuminanceRecord%d.h",LR.screen);
  83.     i=ReadLuminanceRecord(string,&LR,0);        // try to read latest LuminanceRecord
  84.     if(i<1)printf("Warning: couldn't find “%s”. Calibrating screen %d.\n"
  85.         ,string,LR.screen);
  86.     else oldScreen=LR.screen;
  87.  
  88.     mtfStyle[0].continuing=0;
  89.     mtfStyle[0].lineWidth=1;
  90.     mtfStyle[0].symbolWidth=0;
  91.     mtfStyle[0].dash[0]=0;
  92.     mtfStyle[0].dashOffset=0;
  93.     mtfStyle[0].color=blackColor;
  94.     mtfStyle[1]=mtfStyle[0];
  95.     mtfStyle[1].color=blueColor;
  96.     
  97.     GetDateTime(&seconds);
  98.     sprintf(outName,"MTF%d.%s",LR.screen,DatedString(seconds));
  99.     outfile[0]=stdout;
  100.     outfile[1]=fopen(outName,"w");
  101.     if(outfile[1]==NULL) PrintfExit("Sorry, can't create \"%s\".\007\n",outName);
  102.     SetFileInfo(outName,'TEXT','QKPT');    /* for Kaleidagraph */
  103.  
  104.     /* Find device corresponding to the experimental screen. */
  105.     oldGDHandle=GetGDevice();
  106.     device = GetScreenDevice(LR.screen);
  107.     if(NewPaletteManager())
  108.         SetDepth(device,8,1,1);    /* 8-bit pixels, color mode */
  109.     window = GDOpenWindow(device);
  110.     ShowWindow((WindowPtr)window);
  111.     
  112.     if(oldScreen==LR.screen){
  113.         printf("Luminance was calibrated on %s\n",LR.date);
  114.         printf("at %.2f %s by %s\n",LR.LBackground,LR.units,LR.notes);
  115.         printf("Range is %.1f to %.1f %s.\n",LR.LMin,LR.LMax,LR.units);
  116.     }
  117.     L=LR.LBackground;    /* desired mean luminance */
  118.     printf("Enter desired display luminance in cd/m^2 (%.2f):",L);
  119.     fgets(string,sizeof(string)-1,stdin);
  120.     sscanf(string,"%lf", &L);
  121.  
  122.     c=0.9;
  123.     c=MIN(c,MIN((LR.LMax-L)/L,(L-LR.LMin)/L));
  124.     printf("Contrast %.1f\n",c);
  125.     printf("Computing sinusoidal lookup table . . .\n");
  126.     for(i=0;i<256;i++){
  127.         SetLuminance(device,&LR,i,L*(1.0+c*sin(TWO_PI*i/256.0)),L*(1.0-c),L*(1.0+c));
  128.     }
  129.  
  130.     tPeriod=64;                                    /* frames in one temporal period */
  131.     width=window->portRect.right;                /* pixels across screen */
  132.     SetRect(&r,0,0,2*tPeriod,2*tPeriod);
  133.     OffsetRect(&r,64,64);
  134.     plotWindow=NewWindow(NULL,&r,"\pL",1,noGrowDocProc,(WindowPtr) -1L,FALSE,0L);
  135.     SetRect(&r,0,0,192,192);
  136.     OffsetRect(&r,640-r.right-64,64);
  137.     mtfWindow=NewWindow(NULL,&r,"\pMTF",1,noGrowDocProc,(WindowPtr) -1L,FALSE,0L);
  138.     MeasureContrast(device,&LR,tPeriod,cH,plotWindow);        /* measure vBlack */
  139.     ffprintf(outfile,
  140.         "notes" "\tcycles/pixel" "\tcycles/screen" 
  141.         "\tHoriz. Grating Gain" "\tVert. Grating Gain" "\tVert./Horiz.");
  142.     fprintf(outfile[0],"\n");
  143.     fprintf(outfile[1],"\tc" "\tcH[0]" "\tcH[1]" "\tcH[2]" "\tcV[0]" "\tcV[1]" "\tcV[2]" 
  144.         "\t(cV[0]-cH[0])/cH[0]" "\tcV[2]/(c*V/H)^2" "\t(cV[0]-cH[0])/cH[0]/(c*V/H)^2"
  145.         "\n");
  146.     ffprintf(outfile,"notes");    /* put some text in notes column */
  147.     for(fNominal=0.0;;fNominal=MAX(fNominal*pow(2.0,1.0/2.0),fNominal+1.0/width)){
  148.         f=floor(fNominal*width+0.5)/width; /* integral periods for zero mean over line */
  149.         if(fNominal>0.5)f=0.5;
  150.         /* horizontal grating */
  151.         srcRect=dstRect=window->portRect;
  152.         dstRect.left=srcRect.right=1;
  153.         for(i=0;i<srcRect.bottom;i++){
  154.             row[0]=128.0+127.5*saw(TWO_PI*f*i);
  155.             SetWindowPixelsQuickly((WindowPtr)window,0,i,row,1);
  156.         }
  157.         /* CopyBits uses the inverse color table of the current GDevice. */
  158.         SetGDevice(device);
  159.         CopyBits((BitMap *)*window->portPixMap,(BitMap *)*window->portPixMap,
  160.             &srcRect,&dstRect,srcCopy,NULL);
  161.         SetGDevice(oldGDHandle);
  162.         cOut=MeasureContrast(device,&LR,tPeriod,cH,plotWindow);
  163.         ffprintf(outfile,"\t%10.5f" "\t%10.2f" "\t%12.6f",f,f*dstRect.right,cOut/c);
  164.         /* draw the MTF */
  165.         PlotXY(mtfWindow
  166.             ,log(f/(1./640.))/log(0.5/(1./640.)), log(cOut/c/0.3)/log(1.1/0.3)
  167.             ,&mtfStyle[0]);
  168.  
  169.         /* vertical grating */
  170.         srcRect=dstRect=window->portRect;
  171.         dstRect.bottom=srcRect.top=srcRect.bottom-1;
  172.         for(i=0;i<srcRect.right;i++)row[i]=128.0+127.5*saw(TWO_PI*f*i);
  173.         SetWindowPixelsQuickly((WindowPtr)window,0,srcRect.top,row,srcRect.right);
  174.         SetGDevice(device);
  175.         CopyBits((BitMap *)*window->portPixMap,(BitMap *)*window->portPixMap,
  176.             &srcRect,&dstRect,srcCopy,NULL);
  177.         SetGDevice(oldGDHandle);
  178.         cOut=MeasureContrast(device,&LR,tPeriod,cV,plotWindow);
  179.         ffprintf(outfile,"\t%20.6f" "\t%18.6f",cOut/c,cV[1]/cH[1]);
  180.         fprintf(outfile[0],"\n");
  181.         fprintf(outfile[1],"\t%f" "\t%f" "\t%f" "\t%f" "\t%f" "\t%f" "\t%f"
  182.             ,c,cH[0],cH[1],cH[2],cV[0],cV[1],cV[2]);
  183.         a=(cV[0]-cH[0])/cH[0];    /* dc error, expressed as a contrast */
  184.         b=cV[1]/cH[1];            /* contrast gain of video amplifier */
  185.         fprintf(outfile[1],"\t%f" "\t%f" "\t%f" "\n"
  186.             ,a,cV[2]/(c*b*c*b),a/(c*b*c*b));
  187.         /* draw the MTF */
  188.         PlotXY(mtfWindow
  189.             ,log(f/(1./640.))/log(0.5/(1./640.)), log(cOut/c/0.3)/log(1.1/0.3)
  190.             ,&mtfStyle[1]);
  191.         
  192.         if(f==0.5)break;
  193.     }
  194.     
  195.     /* print notes */
  196.     IUDateString(seconds,longDate,(unsigned char *)string);
  197.     p2cstr((unsigned char *)string);
  198.     IUTimeString(seconds,FALSE,(unsigned char *)string2);
  199.     p2cstr((unsigned char *)string2);
  200.     ffprintf(outfile,"%s %s\n",string2,string);
  201.     ffprintf(outfile,"%.1f Hz\n",1.0/(tPeriod*0.015));
  202.     ffprintf(outfile,"%.2f cd/m^2\n",L);
  203.     ffprintf(outfile,"%.3f contrast\n",c);
  204.     ffprintf(outfile,"\n" "screen %d\n",LR.screen);
  205.     ffprintf(outfile,"\"%s\" monitor %s\n",LR.name,LR.id);
  206.     ffprintf(outfile,"on %s\n",LR.date);
  207.     ffprintf(outfile,"by %s\n",LR.notes);
  208.     ffprintf(outfile,"at %.2f %s\n",LR.LBackground,LR.units);
  209.  
  210.     GDDisposeWindow(window);
  211.     if(outfile[1] != NULL)fclose(outfile[1]);
  212.     printf("The sum of the times for GetVoltage() and LoadLuminances() should\n"
  213.         "be about one frame, 15 ms. If it's longer, e.g. 30 ms, you should reduce\n"
  214.         "the value of \"frames\" in MeasureContrast().\n");
  215.     DisposeWindow(plotWindow);
  216.     DisposeWindow(mtfWindow);
  217. }
  218.  
  219.  
  220. double MeasureContrast(GDHandle device,LuminanceRecord *LP
  221.     ,int tPeriod,double c[10],WindowPtr plotWindow)
  222. /*
  223. Measures the contrast of a drifting grating of unknown phase.
  224. The temporal period is tPeriod frames. The c[] array is
  225. filled with the amplitude spectrum, where c[i] is the contrast of the i-th harmonic,
  226. except that c[0], which would always be 1, instead is set to the dc level.
  227. */
  228. {
  229.     double v[256];
  230.     double p[10];
  231.     register int i,j;
  232.     ColorSpec doubleTable[512];    // was static, dgp 6/15/94
  233.     static double sinTable[256],cosTable[256];
  234.     static int tablePeriod=0;
  235.     register double a,b;
  236.     double frequency=2000.,gain=100.;
  237.     double frames=0.6;    /* How long the A/D should spend sampling the photometer */
  238.     long n;
  239.     static double vBlack=0.0;
  240.     static int blackSet=0;
  241.     char string[80];
  242.     WindowPtr oldPort;
  243.     int v0,ip;
  244.  
  245.     assert(StackSpace()>4000);
  246.     if(tPeriod>256){
  247.         printf("Warning. Reducing tPeriod to 256.\n");
  248.         tPeriod=256;
  249.     }
  250.     if(tablePeriod != tPeriod){
  251.         for(i=0;i<tPeriod;i++){
  252.             sinTable[i]=sin(i*TWO_PI/tPeriod);
  253.             cosTable[i]=cos(i*TWO_PI/tPeriod);
  254.         }
  255.         tablePeriod=tPeriod;
  256.     }
  257.     
  258.     n=(long)floor(0.5+frequency*0.015*frames);
  259.  
  260.     RemeasureContrast:
  261.     if(!blackSet){
  262.         doubleTable[0].rgb.red=doubleTable[0].rgb.green=doubleTable[0].rgb.blue=0;
  263.         for(i=0;i<256;i++)doubleTable[i]=doubleTable[i+256]=doubleTable[0];
  264.         LoadLuminances(device
  265.             ,(LuminanceRecord *)(doubleTable+(i%tPeriod)*(256/tPeriod)),0,255);
  266.         printf("Please block all light to set black. Hit cr when ready:");
  267.         gets(string);
  268.         vBlack=0.0;
  269.     }
  270.     else for(i=0;i<256;i++)doubleTable[i]=doubleTable[i+256]=LP->table[i];
  271.  
  272.     /* warm up for one period before collecting data */
  273.     SetPriority(7);
  274.     for(i=0;i<tPeriod;i++){
  275.         LoadLuminances(device
  276.             ,(LuminanceRecord *)(doubleTable+(i%tPeriod)*(256/tPeriod)),0,255);
  277.         GetVoltage(1,&gain,&frequency,n,NULL);
  278.         v[i]=0.0;
  279.     }
  280.     #if SYMANTEC_C_PROFILE
  281.         _profile=1;
  282.     #endif
  283.     for(i=0;i<tPeriod*4;i++){
  284.         LoadLuminances(device
  285.             ,(LuminanceRecord *)(doubleTable+(i%tPeriod)*(256/tPeriod)),0,255);
  286.         v[i%tPeriod]+=GetVoltage(1,&gain,&frequency,n,NULL);
  287.     }
  288.     #if SYMANTEC_C_PROFILE
  289.         _profile=0;
  290.     #endif
  291.     SetPriority(0);
  292.     a=0.0;
  293.     for(i=0;i<tPeriod;i++) a+=v[i];
  294.     c[0]=a/tPeriod-vBlack;
  295.     if(!blackSet){
  296.         vBlack=c[0];
  297.         blackSet=1;
  298.         printf("Black %g mV\n",vBlack*1000.0);
  299.         printf("Now please remove light block. Hit cr when ready:");
  300.         gets(string);
  301.         goto RemeasureContrast;
  302.     }
  303.     for(j=1;j<10;j++){
  304.         a=b=0.0;
  305.         for(i=0;i<tPeriod;i++){
  306.             a+=sinTable[i*j%tPeriod]*v[i];
  307.             b+=cosTable[i*j%tPeriod]*v[i];
  308.         }
  309.         c[j]=2.0*sqrt(a*a+b*b)/tPeriod/c[0];
  310.         p[j]=atan2(a,b);
  311.     }
  312.     /* Show one period of raw data, fundamental, raw minus fundamental, 2nd harmonic */
  313.     GetPort(&oldPort);
  314.     SetPort(plotWindow);
  315.     BringToFront(plotWindow);
  316.     EraseRect(&plotWindow->portRect);
  317.     v0=plotWindow->portRect.bottom/4;
  318.     SetOrigin(0,-2*v0);
  319.     MoveTo(0,-v0*(v[0]-vBlack)/c[0]);
  320.     for(i=1;i<tPeriod;i++)LineTo(i,-v0*(v[i]-vBlack)/c[0]);
  321.     SetOrigin(0,-4*v0);
  322.     j=1;
  323.     ip=tPeriod*(1.0-p[j]/TWO_PI);
  324.     a=-v0*(0.0+c[j]*cosTable[(i*j+ip)%tPeriod]);
  325.     MoveTo(0,-v0*(v[0]-vBlack)/c[0]-a);
  326.     for(i=1;i<tPeriod;i++){
  327.         a=-v0*(0.0+c[j]*cosTable[(i*j+ip)%tPeriod]);
  328.         LineTo(i,-v0*(v[i]-vBlack)/c[0]-a);
  329.     }
  330.     ForeColor(blueColor);
  331.     for(j=1;j<3;j++){
  332.         SetOrigin(-tPeriod,-2*v0*j);
  333.         ip=tPeriod*(1.0-p[j]/TWO_PI);
  334.         i=0;
  335.         a=-v0*(1.0+c[j]*cosTable[(i*j+ip)%tPeriod]);
  336.         MoveTo(0,a);
  337.         for(i=1;i<tPeriod;i++){
  338.             a=-v0*(1.0+c[j]*cosTable[(i*j+ip)%tPeriod]);
  339.             LineTo(i,a);
  340.         }
  341.     }
  342.     ForeColor(blackColor);
  343.     SetOrigin(0,0);
  344.     SetPort(oldPort);
  345.     return c[1];
  346. }
  347.  
  348.  
  349. double saw(double x)
  350. /* returns sawtooth function with same phase, amplitude, and symmetry as sin() */
  351. {
  352.     x/=TWO_PI;
  353.     x-=0.5;
  354.     x-=floor(x);
  355.     return 2.0*x-1.0;
  356. }
  357.